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CONTRACTOR REPORT 


NUMERICAL ANALYSIS OF LAMINAR AND TURBULENT INCOMPRESSIBLE 
FLOWS USING THE FINITE ELEMENT METHOD 

1. INTRODUCTION 


Applications of finite element methods to incompressible viscous flows were 
introduced in the early 1970s [1,2]. In dealing with incompressible flow problems via 
velocity-pressure formulations, there exist two major obstacles to overcome: (1) the 

generation of spurious pressure modes by the coupling of velocity and pressure, and 
(2) oscillatory behaviors (wiggles) of velocity in convection- dominated flows. 

Since the mid-1970s, it was recognized that the mixed interpolation (i.e., the 
interpolation function for pressure as a polynomial of, at least, one order lower than 
that for velocity) is required to prevent the generation of the spurious pressure 
modes [3], Based on the similar motivation, the penalty finite element method was 
introduced in the late 1970s [4,5]. With this method, the pressure term in the momen- 
tum equations could be eliminated by penalizing the incompressibility with a small 
perturbation and, consequently, the number of equations to be solved could be 
reduced. As mentioned in Bercovier and Engelman [4] and Hughes et al. [5], the 
penalty term in momentum equations must be integrated by the reduced integration 
method to prevent overconstraints of incompressibility. The mathematical investigation 
of the penalty finite element method with reduced integration has been performed by 
Oden et al. [6]. As an important improvement of the penalty finite element method, 
Sani et al . [ 7] showed that a concept of the mixed interpolation method could be 
applied to the penalty finite element method. This was done by a consistent applica- 
tion of Galerkin method to the appropriate continuity equation and a consistent inte- 
gration of all terms in momentum equations. This method is called the penalty finite 
element integration with a consistent integration. Engelman et al. [8] showed that 
the reduced integration scheme corresponds to the consistent one only in simple cases, 
and the consistent integration could provide more accurate solutions than the reduced 
integration in most situations. 

When mixed interpolation polynomials for velocity and pressure are used, a 
selection of appropriate polynomials is very important [7]. A comparison of various 
mixed interpolation methods was performed by Huyakorn et al. [9] and they concluded 
that the bilinear interpolations for velocity and piecewise constant pressure, which is 
the simplest combination in two-dimensional domain, generates erroneous pressure 
modes in some cases. Oden and Jacquotte [10] proved both mathematically and 
numerically that the combination of biquadratic velocity with linear pressure is one of 
the most stable choices for smooth pressure modes, and the same conclusion was pro- 
vided by Engelman et al. [8] after numerical experiments in both two- and three- 
dimensional domain. 

It is well-known that the application of conventional Galerkin finite element 
method to convection- dominated flow problems leads to the centered- difference treat- 
ment of the convection term which is a cause of the oscillatory behaviors (wiggles) 
of velocity or any other flow variables. As mentioned in Gresho and Lee [11], the 
obvious way to eliminate wiggles is the mesh refinement in the convection- dominated 
region such that the local Reynolds (or Peclet) number becomes small. On the other 
hand, the wiggle-free solution could also be obtained by the use of the "upwind" 


finite element scheme as an alternative choice of a mesh refinement. The basic idea 
of upwind finite element method was presented by Christie et al. [12] using modified 
weighting functions. In this method, the upstream contribution of the weighting 
function is weighted heavier than that of downstream. This concept was applied to 
2-D cases by Heinrich et al. [13,14], Hughes and Brooks [15] pointed out that the 
above upwinding scheme could generate excessive numerical diffusion perpendicular 
to the flow direction in multi- dimensional cases. To overcome this shortcoming, they 
developed the streamline upwind /Petrov Galerkin method by modifying the weighting 
function to add numerical diffusion only in the flow direction. The same concept of 
streamline upwinding has been used by Gresho et al. [16]. Instead of modifying the 
weight function, they defined an additional term, so-called "balancing tensor diffusiv- 
ity , " which was originated from the second-order difference term in the explicit time 
integration. They also proved that this term influences only in the flow direction in 
both steady and unsteady cases. 

FIDAP is a finite element code for the analysis of incompressible fluid flows and 
heat transfer in a multidimensional domain. This code includes both of the available 
methods for the treatment of velocity- pressure coupling, i.e., penalty and mixed 
interpolation methods with variable choices of interpolation polynomials for velocity 
and pressure. Streamline upwind method is included as an option for the convection- 
dominated flow situations. 


2. TEST CASES 


The finite element method described in the previous section is applied to the 
following four examples, i.e.: (1) 2-D laminar flow inside a wall-driven cavity, (2) 

2-D laminar flow over a backward -facing step, (3) 2-D turbulent flow over a back- 
ward-facing step, and (4) 2-D turbulent flow through a turn-around-duct, are 
selected. The first three examples are well-known benchmark problems which are 
useful for the evaluation of any CFD methodologies and there exists computational 
and/or experimental data for these problems. The fourth case is chosen as an appli- 
cation of the present method to the flow inside a Space Shuttle Main Engine (SSME). 

For all test cases, the nine-nodes quadrilateral finite element with linear dis- 
continuous interpolations for pressure (Q 2 /P 1 ) is used. This approach is known to 

be one of the optimal choices for the treatment of velocity-pressure coupling of the 

- fi 

incompressible fluid flow. And, the selected value of the penalty parameter e is 10 
Also, the influence of the streamline upwinding (STU) method and its necessary has 
been investigated. 

2-1. 2-D Laminar Flow Inside a Wall-Driven Cavity 

As the first test case, a 2-D laminar flow inside a wall-driven cavity is chosen. 
The essential point of this test case is the prediction of various vortices inside a 
cavity as shown in Figure 1. The computational solutions of the present study are 
compared to those worked by Ghia et al. [18] who used the finite difference method 
of streamline- vorticity formulations with fine (129 x 129 grid points for Re < 3200 
and 257 x 257 grid points for Re > 5000) but uniform meshes. 

In the present study, 40 x 40 elements (81 x 81 grid points) are used. Fine 
meshes are used near walls because of viscous boundary layers and the secondary/ 
tertiary recirculating zones of small sizes near corners. The characteristic length of 
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the smallest element is 0.00326 at four corners and the largest one is 0.03074 at the 
center of a cavity. The selected Raynolds numbers are 100, 400, 1000, 3200, 5000, 

7500, and 10000. To minimize CPU time and optimize initial conditions of the high 
Reynolds number flow, the restart procedure with the increment of Reynolds number 
is used. In case of Re = 100, the solution of Stokes flow is chosen as the initial 
condition. At every selected Reynolds number, the first three iterations for nonlinear 
solutions are done by the successive substitution method and the quasi- Newton method 
is chosen after fourth iteration. With these combinations, solutions converged smoothly 
to 1 percent convergence criteria of the relative velocity, within 4 to 5 iterations at 
each step. When STU was selected, the number of iterations were unchanged for all 
Reynolds numbers, but the slope of the convergence curve was steeper than the other. 

Figure 2 shows resulting streamline contours for various Reynolds numbers. 
Influence of STU is also shown at each Reynolds number. Labels of streamline con- 
tours are shown in Table 1. Results without STU are generally comparable to those 
of Ghia et al. with locations and shapes of secondary and tertiary vortices near the 
bottom corners. As expected, STU generates strong numerical diffusions particularly 
in high Reynolds number flows. These numerical diffusivities reduce sizes of various 
vortices. Particularly, the primary vortex in the center of a cavity is seriously 
influenced by STU. 

Pressure contours are shown in Figure 3 with labels in Table 1. STU reduces 

total pressure differences (AP = P -P . ) in all Reynolds numbers. Comparing 

max min ° 

those results with the results of Gresho et al. [19] for Re = 5000 and 10000, the AP's 

of the present study without STU are bigger than those obtained by Gresho et al., 

but similar trends of pressure contours are generated in both cases. 

Horizontal velocities along the vertical centerline inside a cavity are compared 
with Ghia et al. in Figure 4. In cases without STU, velocity profiles in all Reynolds 
numbers excellently agree with Ghia et al. The minimum difference of the peak 
velocity is 1.00 percent at Re = 100 and the maximum is 1.57 percent at Re = 10000. 

On the other hand, when STU is selected, the minimum difference is 8.241 percent 
at Re = 100 and the maximum is 26.10 percent at Re = 10000. Overall velocity pro- 
files with STU are quite different from those without STU at very high Reynolds 
numbers . 

Extreme values of stream functions at various vortices inside a cavity are shown 
in Table 2. Results of the present study (both with and without STU) are consis- 
tently lower than those of Ghia et al. in all Reynolds numbers. For the cases without 
STU, differences from Ghia et al. are 1.03 percent to 6.40 percent for the primary 
vortex, 3.28 percent to 18.1 percent for the bottom vortices, and 9.85 percent to 
16.67 percent for the top left vortex. Values of extreme stream functions of the 
primary vortex are higher than those of Gresho et al. This may be simply due to 
finer grids (h = 0.03074) employed in the central region than Gresho et al. (h = 0.060). 
Another difference between two cases is the underprediction of bottom vortices of the 
present study and the overprediction of Gresho et al. Once again, it is shown that 
STU seriously influences the strengths of vortices. 


2-2. 2-D Laminar Flow Over a Backward-Facing Step 


Two-dimensional laminar flow over a backward- facing step is selected as the 
second example. Geometry and boundary conditions are shown in Figure 5. The 
aspect ratio of the backward- facing step (h) to the overall cross-sectional width is 
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1:2 and the total length in the horizontal direction is 30 h. A fully- developed para- 
bolic velocity profile is prescribed at the inlet boundary. Predictions of recirculating 
zones behind a backward- facing step and near the upper wall are compared with 
experimental measurements of Armaly et al. [20]. 

30 x 16 elements (61 x 33 grid points) are selected with fine meshes near walls 
and mixing zone along the centerline. Selected Reynolds numbers are 100, 200, 300, 
400, 500, 600, 700, and 800. Here, the Reynolds number is defined by the bulk 
velocity at the inlet boundary and cross-sectional width of the whole domain which is 
the same definition as Armaly et al. The same restart procedure as the previous 
case is adopted. Also, the combination of the successive substitution and quasi- 
Newton methods is used to stabilize iterations for the nonlinear solutions. Four to 
five iterations were performed in each step for 1 percent convergence of relative 
velocity and the same trends of convergence as the previous case were observed. 

Figure 6 represents streamline contours for selected Raynolds numbers. Labels 
of contours are shown in Table 3. The recirculating zone behind a backward- facing 
step becomes bigger as the Reynolds number increased, and another recirculating 
zone near the upper wall is generated for Re = 500. STU reduces streamline levels 
of both recirculating zones. But the reattachment length of the recirculating zone 
behind a backward- facing zone, x^, is a little longer than the other case for all 
Reynolds numbers. 

Pressure contours are shown in Figure 7 and Table 3 for selected Reynolds 
numbers. As Reynolds number increases, the order of magnitude of pressure becomes 
smaller and the position of maximum pressure moves to downstream. At a fixed 
Reynolds number, the total pressure difference in the case with STU is bigger than 
that without STU. 

Comparisons with experimental data of Armaly et al. are shown in Figure 8. 
There is excellent agreement between predicted reattachment lengths and experimental 
measurements. Discrepancies between computational predictions and experimental data 
at high Reynolds numbers are originated by the 3-D effects of experiments. The 
general trend of the curve of x^ (the reattachment length of the recirculating zone 

behind a backward-facing step) is the same as that of Kim and Moin [21] who used 
the finite difference method in the calculations with 101 x 101 grid points. Also, the 
prediction of beginning (x 2 ) and ending (x^) points of the recirculating zone near 

the upper wall compares satisfactorily with experimental data. By using STU, x^ 

becomes longer from 1.3 percent (Re = 800) to 8.0 percent (Re = 100) and x 2 and x^ 

are shifted to downstream at low Reynolds number and upstream at high Reynolds 
number. The width of the recirculating zone near the upper wall is decreased from 
a maximum of 34.4 percent at Re = 500 and to a minimum of 2.7 percent at Re = 800 
by STU. 


2-$. 2-D Turbulent Flow Ove** a Backward-Facing Step 

A two-dimensional turbulent flow over a backward- facing step, which is one of 
the standard test cases of complex turbulent flows presented in the Stanford Confer- 
ence [22], is chosen. The geometry and boundary conditions are shown in Figure 9. 
The aspect ratio of the backward- facing step (h) to the whole cross-sectional width 
is 1:3. Total length in the horizontal direction behind the backward- facing step is 
24 h. As previously mentioned, the standard k- e turbulence model of Launder and 
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Spalding [17] is used with wall- functions at non-slip walls. Uniform velocity is given 
at the inlet boundary. Inlet boundary conditions for k and e are prescribed as 
follows: 


k = 0.003 u Q 2 


e 



.1.5 


1 


-1 


where 1 = 0.03 h. 

The Reynolds number based on the inlet velocity and the height of a backward- 
facing step is chosen as 69,610 which represents the same flow condition as Kim et al. 
[23]. 


22 x 16 nine-nodes quadrilateral elements (45 x 33 grid points) are used with 
fine meshes near walls. For the iterative solutions of nonlinear equations, the succes- 
sive substitution method (acceleration factor = 0.5) is selected. When STU is used, 
solution was converged to the convergence criteria 0.001 percent of the relative 
velocity with 65 iterations. Newton-type iterative methods were not suitable to this 
particular problem because of complex nonlinear terms in the turbulence equations. 

The streamline upwind (STU) scheme was necessary. When STU was not selected, 
wiggles are generated in the solution of the e- equation after 4 to 5 iterations, and 
velocities are distorted by the influence of these wiggles. 

Figure 10 represents streamline contours. The most important parameter to 
compare with experimental data [22,23] is the reattachment length, x f , of the separa- 
tion zone behind a backward-facing step. Predicted x p in the present study is 5.59. 
The acceptable value of x p in the experimental measurement is 7.0 ± 1.0 [23]. 
Generally, k-e turbulence model underpredicts x r by about 20 percent as mentioned 

in Eaton [24]. In most computational results with k-e turbulence model, the predicted 
is in the range of 5.2 to 5.8 [ 25 ]. The exact value of x^. depends on the number 

of grid points and computational methodologies. Also, it depends on various turbu- 
lence models [ 24] . Predicted x^ in the present study is comparable to other computa- 
tional results with k-e turbulence model. In addition, although it is not clearly 
shown in Figure 10, a very small but apparent secondary recirculation exists near 
the bottom corner of a backward- facing step. This was also predicted by Kaul [26]. 

In Figure 11, the predicted mean velocity profiles are compared with experi- 
mental measurements. Predictions reasonably agree with experimental data downstream, 
but there exist some differences between predictions and measurements in the separa- 
tion region and near the wall. Figure 12 shows a comparison of pressure coefficients 
C . In the present study, velocity and pressure at the center of inlet boundary are 

hr 

selected as reference velocity and pressure for the calculation of pressure coefficients. 
Cp is generally predicted quite well except in the separation zone. The predicted 

value of Cp in downstream is 0.36. 
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Contours of turbulent kinetic energy (TKE) are shown in Figure 13. The peak 
TKE exist at the center of a separation region and steep gradients of TKE exist near 
the centerline and solid walls. Labels of streamline and TKE contours are shown in 
Table 4. Predicted TKE's are compared with experimental data in Figure 14. The 
location of maximum TKE moves closer to the wall in downstream. Differences between 
predictions and measurements are apparent near the separation and mixing zones. 

The overprediction of TKE in the separation region may be the cause of the under- 
prediction of the reattachment length. Similar trends of the overprediction of TKE 
are experienced in other computational results using the k- e turbulence model [24,25]. 


2-4. 2-D Turbulent Flow Through a Turn-Around-Duct 

This example is selected as an application of the present study to the design of 
flow paths inside a Space Shuttle Main Engine (SSME). Transition duct between the 
fuel-side preburner and hot-gas manifolds inside SSME is a 180-degree curved annular 
turn- around- duct. Experimental studies for turbulent flow through this geometry are 
done by Sandborn [27] and Sharma and Ostermier [28] using rectangular and axi- 
symmetric shape turn-around- ducts. Experimental apparatus for both cases are shown 
in Figure 15. Computational results of the present study are compared with these 
experimental cases. 

In case of a rectangular turn- around- duct , the width of a duct (h) is 1.0 and 
the radius of curvature along the centerline of a duct is 1.0. As shown in Figure 16, 
two kinds of computational domains are selected. As a first case, the inlet boundary 
is set at the same location as Sandborn and the uniform velocity (U| n = 1*0) has been 

prescribed as the inlet boundary condition. As the other case, the inlet boundary is 
set at 5.2 h ahead of the turn and the measured data of mean velocities of Sandborn 
are prescribed as the inlet boundary condition. The purpose of this study is to 
investigate the influence of inlet boundary conditions. The wall function is applied 
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along solid boundaries in both cases. And, Reynolds number is selected as 10 . 

40 x 14 finite elements (81 x 29 grid points) and 35 x 14 finite elements (71 x 29 grid 
points) are used in cases 1 and 2, respectively. The streamline upwind (STU) was 
necessarily selected to stabilize numerical solutions. 

Computational results of cases 1 and 2 are shown in Figures 17 and 18, respec- 
tively. Velocity profiles near the inner wall after a turn in case 2 are almost sepa- 
rated from the wall, but no recirculating zone is predicted in both cases which was 
observed in Sandborn’s experiments. This may be originated by the contribution of 
strong numerical diffusions due to STU in momentum equations. Pressure gradients 
along the inner wall are steeper than the outer wall in both cases. Gradients of 
TKEs near the inner wall are steeper than the outer wall in both cases but case 2 has 
steeper gradients than case 1. Predicted mean velocity profiles of both cases are 
compared with experimental measurements at various locations in Figure 19, and also, 
with Chen and Sandborn [29] who used a finite difference method in the same geo- 
metry with the second order upwind for momentum equations and a skew upwind for 
k-e equations. Difference between predictions and measurements become greater 
downstream. As expected, predictions of case 2 are better than those of case 1. 
Pressure coefficients along inner and outer walls are compared with experimental 
measurements in Figure 20. Both cases are reasonably comparable with experiments. 

Computational domain of an axisymmetric shape turn- around- duct is shown in 
Figure 21. 40 x 14 finite elements (81 x 29 grid points) are used with STU. Uni- 

form velocity (U_ n = 1.0) is prescribed at the inlet boundary and wall-function is 
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used along the solid wall. Computational results are shown in Figure 22. Velocity 
profiles near the inner wall after a turn are physically more unstable than those of 
a rectangular turn- around- duct . The peak point of the pressure along the inner 
wall is moved to upstream and gradients of TKEs near the inner wall after a turn are 
steeper than those of a rectangular turn-around-duct. Predicted pressure coefficients 
of the present study are compared with experimental measurements of Sharma and 
Ostermier [28] and computational results of Chen [30] which used higher-order 
accurate up winding for momentum equations. Results of the present study are reason- 
ably comparable with experimental data in the upstream of the turn. Differences 
between predictions of the present study and experimental data in the downstream 
region may be originated by two main reasons, i.e., (1) numerical diffusions in 
momentum equations and (2) non-existence of the curvature effects in k-e turbulence 
model. As shown in the results of Chen, predictions can be improved by using 
higher order upwind in momentum equations and modifying k- e turbulence model to 
include curvature effects (extended k-e turbulence model). 


3. CONCLUSIONS 


Some classical benchmark problens in 2-D laminar and turbulent incompressible 
fluid flows are tested for the evaluation of numerical accuracy of the finite element, 
FIDAP , method. As a result of the present study, the following conclusions are drawn: 

1. The consistent integration penalty finite element method with elements 

is efficient to treat velocity-pressure couplings of the incompressible fluid flow. 

2. Streamline upwinding (STU) produces too strong numerical diffusions in the 
high Reynolds number laminar flow, particularly, inside a wall- driven cavity. 

3. When STU is not selected, numerical results of the present study on laminar 
flow are in good agreement with other CFD solutions and experimental data. 

4. STU is essential to the computation of turbulent flow because it stabilizes 
numerical solutions of k-e equations. However, as shown in the case of the flow 

through a turn- around- duct , momentum equations are still strongly influenced by STU. 
Therefore, numerical diffusions in the velocity field must be minimized by applying 
higher-order upwinding techniques to momentum equations. 

5. Results of the present study on turbulent flow correspond to the general 
trends of other computational solutions of the standard k-e turbulence model. Dif- 
ferences from experimental data in both test cases for turbulent flow are mainly 
caused by the limitations of the standard k- e turbulence model on the recirculating 
and curved flow. Therefore, the advanced turbulence models are necessarily required 
for the accurate predictions of complex turbulent flows. The incorporation of such 
advanced turbulence models in the finite element code, FIDAP, is one of several tasks 
requiring a continuous study. 
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Figure 1. Geometry and boundary conditions of 2-D laminar 
flow inside a wall-driven cavity. 



























(WITHOUT STU) 


(WITH STU) 


(a) Re = 100 


Figure 3. Pressure contours of 2-D laminar flow. 
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(WITH STU) 


(e) Re = 5000 


Figure 3. (Continued) 











Re = 10000 Re = 7500 Re = 5000 Re = 3200 Re = 1000 Re = 400 Re =100 


Figure 4. A comparison of u- velocities along a vertical centerline, — : present 

study without STU, : present study with STU, 

O: Ghia et al. [18] . 





Figure 5. Geometry and boundary conditions of 2-D laminar flow 
over a backward- facing step. 
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Figure 7. Pressure contours of 2-D laminar flow over a backward- facing step. 
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Figure 7. (Concluded) 
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Figure 11. A comparison of predicted and measured velocities for 2-D turbulent flow 

over a backward- facing step, : present study, O: experimental 

data of Kim et al. [23] . 




o 


o 


o 


o 



oefficient 


for 








END VIEW OF TURNAROUND BASE 

(b) AXI-SYMMETRIC TURN-AROUND DUCT 


Figure 15. Experimental apparatus of rectangular and 
axisymmetric turn-around ducts. 
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(a) CASE 1 



(b) CASE 2 


Figure 16. Schematics of the rectangular turn-around duct. 










(a) 1.52 AHEAD OF THE TURN 


(b) START OF THE TURN 


Figure 19. 
of 2- 


A comparison of predicted and measured mean velocity profiles 
turbulent flow through a rectangular turn-around duct. 
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Figure 20. A comparison of predicted and measured pressure coefficients of 2-D 
turbulent flow through a rectangular turn-around duct. 
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Figure 23. A comparison of predicted and measured pressure coefficients of 2-D turbulent 

flow through an axisymmetric turn-around duct. 
















TABLE 2. EXTREME VALUES OF STREAM FUNCTIONS AT VARIOUS VORTICES 

INSIDE A WALL-DRIVEN CAVITY 
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TABLE 3. VALUES OF 1 

IN 


LABEL STREAMLINE 

a -0.0250 

b -0.0150 

c -0.0050 

d 0.0000 

e 0.1000 

f 0.2000 

9 0.3000 


0.4000 

0.4391 

0.4408 
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PRESSURE 


Re = 300 

-0.0075 

-0.0050 

-0.0025 

- 0.0000 

-0.0050 

- 0.0100 

-0.0150 

- 0.0200 
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-0.0300 

0.0350 
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Re = 500 

-0.0350 
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0.0075 
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TABLE 4. VALUES OF STREAMLINE AND TURBULENT KINETIC 
ENERGY CONTOURS IN FIGURES 10 AND 11 


LABEL 

STREAMLINE 

TURBULENT 
KINETIC ENERGY 

a 

-0.0600 

0.0450 
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-0.0400 

0.0400 

c 

-0.0200 

0.0350 

d 

0.0000 

0.0300 
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0.2000 

0.0250 
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0.4000 

0.0200 

g 

0.6000 

0.0150 

h 

0.8000 

0.0100 

i 

1.0000 


i 

1 .2000 


k 

1.4000 


i 

1.6000 


m 

1.8000 
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